Case Study: Factors Influencing Life Expectancy using Linear Regression¶
Context:¶
- There have been lot of studies undertaken in the past, on factors affecting life expectancy, considering demographic variables, income composition and mortality rates.
- It was found that affect of immunization and human development index was not taken into account in the past.
- Important immunization like Hepatitis B, Polio and Diphtheria will also be considered. In a nutshell, this case study will focus on immunization factors, mortality factors, economic factors, social factors and other health related factors as well.
- In this case study, we will use linear regression to see the effect of various factors on Life Expectancy.
Problem:¶
The data-set aims to answer the following key questions:
- Does various predicting factors really affect the Life expectancy?
- What are the predicting variables actually affecting the life expectancy?
- Should a country having a lower life expectancy value increase its healthcare expenditure in order to improve its average lifespan?
- Do Infant and Adult mortality rates affect life expectancy?
- Does Life Expectancy has positive or negative correlation with a country's status (developing or developed), lifestyle, GDP, etc.
- What is the impact of schooling on the lifespan of humans?
- Does Life Expectancy have positive or negative relationship with drinking alcohol?
- What is the impact of Immunization coverage (for various disease like Measles,Hepatitis B) on life Expectancy?
Attribute Information:¶
- Country: Country
- Year: Year
- Status: Developed or Developing status
- Life expectancy: Life Expectancy in age
- Adult Mortality: Adult Mortality Rates of both sexes (probability of dying between 15 and 60 years per 1000 population)
- infant deaths: Number of Infant Deaths per 1000 population
- Alcohol: Alcohol, recorded per capita (15+) consumption (in litres of pure alcohol)
- percentage expenditure: Expenditure on health as a percentage of Gross Domestic Product per capita(%)
- Hepatitis B: Hepatitis B (HepB) immunization coverage among 1-year-olds (%)
- Measles: Measles - number of reported cases per 1000 population
- BMI: Average Body Mass Index of entire population
- under-five deaths: Number of under-five deaths per 1000 population
- Polio: Polio (Pol3) immunization coverage among 1-year-olds (%)
- Total expenditure: General government expenditure on health as a percentage of total government expenditure (%)
- Diphtheria: Diphtheria tetanus toxoid and pertussis (DTP3) immunization coverage among 1-year-olds (%)
- HIV/AIDS: Deaths per 1 000 live births HIV/AIDS (0-4 years)
- GDP: Gross Domestic Product per capita (in USD)
- Population: Population of the country
- thinness 1-19 years: Prevalence of thinness among children and adolescents for Age 10 to 19 (% )
- thinness 5-9 years: Prevalence of thinness among children for Age 5 to 9(%)
- Income composition of resources: Human Development Index in terms of income composition of resources (index ranging from 0 to 1)
- Schooling: Number of years of Schooling(years)
Import libraries¶
In [1]:
# Import necessary libraries.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# To enable plotting graphs in Jupyter notebook
%matplotlib inline
Load and explore the data¶
In [2]:
# Load the data into pandas dataframe
data = pd.read_csv('Life Expectancy Data.csv') # Make changes to the path depending on where your data file is stored.
In [3]:
data.head()
Out[3]:
In [4]:
# Check number of rows and columns
data.shape
Out[4]:
In [5]:
# Have a look at the column names
data.columns
Out[5]:
In [6]:
# Check column types and missing values
data.info()
In [7]:
# remove the rows of data which have missing value(s)
data = data.dropna()
In [8]:
# Check the unique values in each column of the dataframe.
data.nunique()
Out[8]:
Insights:¶
- The "Status" column has 2 unique values. i.e. The values are "Developing" and "Developed"
- The "Country" column has 133 unique values. i.e. The data is collected from 133 countries.
In [9]:
plt.figure(figsize=(10,7))
plt.scatter(data['Schooling'], data['Life expectancy'], color='red')
plt.title('Life expectancy Vs Schooling', fontsize=14)
plt.xlabel('Schooling', fontsize=14)
plt.ylabel('Life expectancy', fontsize=14)
plt.grid(True)
plt.show()
Insights:¶
- Clearly, we can see in the above plot that, when the number of years of schooling is increasing, the Life expectancy is also increasing.
- The relationship is linear, as can be seen from the graph.
Dependent vs. Target¶
In [10]:
plt.figure(figsize=(10,7))
plt.scatter(data['Measles'], data['Life expectancy'], color='red')
plt.title('Life expectancy Vs Measles ', fontsize=14)
plt.xlabel('Measles', fontsize=14)
plt.ylabel('Life expectancy', fontsize=14)
plt.grid(True)
plt.show()
In [11]:
plt.figure(figsize=(10,7))
plt.scatter(data['Alcohol'], data['Life expectancy'], color='red')
plt.title('Life expectancy Vs Alcohol', fontsize=14)
plt.xlabel('Alcohol', fontsize=14)
plt.ylabel('Life expectancy', fontsize=14)
plt.grid(True)
plt.show()
In [12]:
plt.figure(figsize=(10,7))
plt.scatter(data['Total expenditure'], data['Life expectancy'], color='red')
plt.title('Life expectancy Vs Total expenditure', fontsize=14)
plt.xlabel('Total expenditure', fontsize=14)
plt.ylabel('Life expectancy', fontsize=14)
plt.grid(True)
plt.show()
In [13]:
plt.figure(figsize=(10,7))
plt.scatter(data['Adult Mortality'], data['Life expectancy'], color='red')
plt.title('Life expectancy Vs Adult Mortality', fontsize=14)
plt.xlabel('Adult Mortality', fontsize=14)
plt.ylabel('Life expectancy', fontsize=14)
plt.grid(True)
plt.show()
Insights:¶
- Here, we see that, as the Adult Mortality is increasing, the Life expectancy is decreasing.
- So, there is inverse relationship between this feature on X-axis and the Target variable on Y-axis.
- Even though the relationship is inverse, but it is linear relationship, as can be seen from the graph.
Conclusion:¶
- Yes, by above analysis, we can conclude that, the dataset is a good fit for the application of Linear Regression.
In [14]:
plt.figure(figsize=(10,7))
plt.scatter(data['Life expectancy'][:200], data['Country'][:200], color='red')
plt.title('Country Vs Life expectancy', fontsize=14)
plt.xlabel('Life expectancy', fontsize=14)
plt.ylabel('Country', fontsize=14)
plt.grid(True)
plt.show()
Insights:¶
- Life expectancy by the country is having a spread of data. (i.e. Some countries might have higher life expectancy and some lower, in the graph, we have plotted 200 data points which are from 15 countries (The names of countries are on Y-axis)
- However, we can see that the red dots are more to the right (near 90 years) for some countries and to the left (less than 50 years) for some countries.
In [15]:
sns.pairplot(data, height=3, diag_kind='auto', corner=True)
plt.show()
Insights:¶
- We can see from the 2nd and 3rd column of the graph that, these features have linear (or inversely linear) relationship.
In [16]:
plt.figure(figsize=(10,10))
sns.boxplot(data['Life expectancy'], orient='v')
plt.show()
Insights:¶
- The bottom black horizontal line of blue box plot is minimum value, which is around 48 years.
- First black horizontal line of rectangle shape of blue box plot is First quartile or 25 percentile. The age corresponding to this line is 64 years.
- Second black horizontal line of rectangle shape of blue box plot is Second quartile or 50% or median. The age corresponding to this is 72 years.
- Third black horizontal line of rectangle shape of blue box plot is third quartile or 75%. The age corresponding to this is 75 years.
- Top black horizontal line of rectangle shape of blue box plot is maximum value. The value corresponding to this is 88 years.
- Small diamond shape of blue box plot is outlier data or erroneous data. Which is below the bottom black horizontal line.
In [17]:
plt.figure(figsize=(8,8))
sns.boxplot(x="Status",y="Life expectancy",data=data)
plt.show()
Insights:¶
- The life expectancy in the developed countries is much higher than the developing countries.
- The median value of life expectancy (approximately as can be seen from the boxplot) for:
- Developing: 69 years
- Developed: 78 years
In [18]:
data[data.columns[:]].corr()['Life expectancy'][:]
Out[18]:
In [19]:
plt.figure(figsize=(20,20))
sns.heatmap(data.corr(), annot=True, fmt=".2")
plt.show()
In [20]:
X = data.drop('Life expectancy', axis=1)
y = data[['Life expectancy']]
print(X.head())
print(y.head())
In [21]:
print(X.shape)
print(y.shape)
In [22]:
X.loc[:, 'Country'] = X.loc[:, 'Country'].astype('category')
X.loc[:, 'Status'] = X.loc[:, 'Status'].astype('category')
X.loc[:, 'Country'] = X.loc[:, 'Country'].cat.codes
X.loc[:, 'Status'] = X.loc[:, 'Status'].cat.codes
X.head()
Out[22]:
In [23]:
X = X.values
y = y.values
In [24]:
#split the data into train and test
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
In [25]:
from sklearn.linear_model import LinearRegression
linearregression = LinearRegression()
linearregression.fit(X_train, y_train)
print("Intercept of the linear equation:", linearregression.intercept_)
print("\nCOefficients of the equation are:", linearregression.coef_)
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
pred = linearregression.predict(X_test)
In [26]:
# Mean Absolute Error
mean_absolute_error(y_test, pred)
Out[26]:
The mean absolute error (MAE) is the simplest regression error metric to understand. We’ll calculate the residual for every data point, taking only the absolute value of each so that negative and positive residuals do not cancel out. We then take the average of all these residuals. Effectively, MAE describes the typical magnitude of the residuals.
In [27]:
# RMSE
mean_squared_error(y_test, pred)**0.5
Out[27]:
The root mean square error (RMSE) is just like the MAE, but squares the difference before summing them all instead of using the absolute value. And then takes the square root of the value.
In [28]:
# R2 Squared:
r2_score(y_test, pred)
Out[28]:
R^2 (coefficient of determination) regression score function.
Best possible score is 1.0 and it can be negative (because the model can be arbitrarily worse). A constant model that always predicts the expected value of y, disregarding the input features, would get a R^2 score of 0.0.
In [29]:
# Training Score
linearregression.score(X_train, y_train)
Out[29]:
In [30]:
# Testing score
linearregression.score(X_test, y_test)
Out[30]:
X_test Vs. Predicted values:¶
In [31]:
df = pd.DataFrame({'Actual': y_test.flatten(), 'Predicted': pred.flatten()})
df
Out[31]:
We can also visualize comparison result as a bar graph using the below script :
Note: As the number of records is huge, for representation purpose I’m taking just 25 records.
In [32]:
df1 = df.head(25)
df1.plot(kind='bar',figsize=(16,10))
plt.grid(which='major', linestyle='-', linewidth='0.5', color='green')
plt.grid(which='minor', linestyle=':', linewidth='0.5', color='black')
plt.show()
- We can observe here that our model has returned pretty good prediction results.
- The Actual and Predicted values are comparable.
LinearRegression model:¶
In [33]:
linearregression_ = LinearRegression(normalize=True) # Linear Regression with hyperparameter normalize=True.
linearregression_.fit(X_train, y_train)
print("Intercept of the linear equation:", linearregression_.intercept_)
print("\nCOefficients of the equation are:", linearregression_.coef_)
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
pred = linearregression_.predict(X_test)
In [34]:
# Mean Absolute Error
mean_absolute_error(y_test, pred)
Out[34]:
In [35]:
# RMSE
mean_squared_error(y_test, pred)**0.5
Out[35]:
In [36]:
# R2 Squared:
r2_score(y_test, pred)
Out[36]:
In [37]:
# Training Score
linearregression_.score(X_train, y_train)
Out[37]:
In [38]:
# Testing score
linearregression_.score(X_test, y_test)
Out[38]:
The Training and testing scores are around 83% and both scores are comparable, hence the model is a good fit.
R2_score is 0.83, that explains 83 % of total variation in the dataset. So, overall the model is very satisfactory.
Setting
normalize=Truedo have impact on co-efficients but they dont affect the best fit line anyway. (refer the coefficients and intercepts are changed, but the accuracy is same as before.)
In [39]:
import statsmodels.api as sm
X = sm.add_constant(X)
linearmodel = sm.OLS(y, X).fit()
predictions = linearmodel.predict(X)
print_model = linearmodel.summary()
print(print_model)
Interpreting the Regression Results:¶
- Adjusted. R-squared:
- The value for Adj. R-squared is 0.837, which is good!
- const coefficient .
- It means that if all the dependent variables (features: like Country, status, Adult mortality and so on..) coefficients are zero, then the expected output (i.e., the Y) would be equal to the const coefficient.
- In our case, the value for const coeff is 308.5849
- Status coeff: It represents the change in the output Y due to a change in the status (everything else held constant).
Schooling coeff: It represents the change in the output Y due to a change of one unit in the Schooling (everything else held constant).
P >|t|:
- A p-value of less than 0.05 is considered to be statistically significant.
- Confidence Interval: It represents the range in which our coefficients are likely to fall (with a likelihood of 95%).
In [40]:
# Plot between residual(actual - predicted) and predicted values
plt.figure(figsize=(10,8))
plt.scatter(linearmodel.resid, linearmodel.predict(), marker='*')
plt.show()
In [41]:
# error distribution
plt.figure(figsize=(10,8))
sns.distplot(linearmodel.resid, hist=True, kde=False, color='red')
plt.show()
